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Numerical relativity has traditionally been pursued via finite differencing. Here we explore pseu- 
dospectral collocation (PSC) as an alternative to finite differencing, focusing particularly on the 
solution of the Hamiltonian constraint — an elliptic partial differential equation — for a black hole 
spacetime with angular momentum and for a black hole spacetime superposed with gravitational 
radiation. In PSC, an approximate solution, generally expressed as a sum over a set of orthogonal 
basis functions (e.g., Chebyshev polynomials), is substituted into the exact system of equations and 
the residual minimized. For systems with analytic solutions the approximate solutions converge 
upon the exact solution exponentially as the number of basis functions is increased. Consequently, 
PSC has a high computational efficiency: for solutions of even modest accuracy we find that PSC is 
substantially more efficient, as measured by either execution time or memory required, than finite 
differencing; furthermore, these savings increase rapidly with increasing accuracy. 
For example, investigating the Hamiltonian constraint equation for a black hole with angular momen- 
tum we find that, where finite difference solutions require a resolution of 1024 x 384 (radial x angular) 
grid points to find a solution of fractional error 10"^ in the ADM mass, a PSC solution achieves 
the same accuracy with only 12 x 4 collocation points. Furthermore, the fractional error is reduced 
to 10"^" by increasing the PSC resolution to 24 x 8, while the same increase in the finite differ- 
ence solution would require (if it were possible) an increase in resolution by a factor of 300 in each 
dimension. Commensurate with the doubling of the resolution in each of the two dimensions the 
computing time required to find the spectral solution increase by a factor of 2^ while the computing 
time required by the finite difference method would increase by a factor of 300^^, or 10^. 
The solution provided by PSC is an analytic function given everywhere, not just at the collocation 
points. Consequently, no interpolation operators need to be defined to determine the function 
values at intermediate points and no special arrangements need to be made to evaluate the solution 
or its derivatives on the boundaries. Since the practice of numerical relativity by finite differencing 
has been, and continues to be, hampered by both high computational resource demands and the 
difficulty of formulating acceptable finite difference alternatives to the analytic boundary conditions, 
PSC should be further pursued as an alternative way of formulating the computational problem of 
finding numerical solutions to the field equations of general relativity. 



I. INTRODUCTION AND SUMMARY 

The partial differential equations (PDE) of numeri- 
cal relativity have typically been solved using finite dif- 
ference methods. In finite differencing (FD) one first 
chooses a finite number of coordinate "grid" points a;„ 
and approximates the space and time derivatives in the 
PDEs by ratios of differences between field and coordi- 
nate values on the grid. With a choice of grid and "dif- 
ferencing scheme" for converting derivatives to ratios of 
differences, the equations of general relativity are approx- 
imated by a system of algebraic equations whose solution 
approximates that of the underlying PDEs. 

In this paper we explore an alternative method for solv- 
ing the elliptic PDEs encountered in numerical relativity: 
pscudospectral collocation (PSC). In PSC one begins by 
postulating an approximate solution, generally as a sum 
over some finite basis of polynomials or trigonometric 



functions. The coefficients in the sum are determined 
by requiring that the residual error, obtained by substi- 
tuting the approximate solution into the exact PDEs, is 
minimized in some suitable sense. Thus, if one describes 
FD as finding the exact solution to an approximate sys- 
tem of equations, one can describe PSC as finding an 
approximate solution to the exact equations. 

Pseudospectral collocation has been applied success- 
fully to solve problems in many fields, including fluid 
dynamics, meteorology, seismology, and relativistic as- 
trophysics (cf. 0-^). Its advantage over FD arises for 
problems with smooth solutions, where the approximate 
solution obtained using PSC converges on the actual so- 
lution exponentially as the number of basis functions is 
increased. The approximate FD solution, on the other 
hand, never converges faster than algebraically with the 
number of grid points. While the computational cost 
per "degree of freedom" — basis functions for PSC, grid 
points for FD — is higher for PSC than for FD, the com- 
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putational cost of a high accuracy PSC solution is a small 
fraction of the cost of an equivalent FD solution. Even for 
problems in which only modest accuracy is needed, PSC 
generally results in a significant computational savings in 
both memory and time compared to FD, especially for 
multidimensional problems. 

The elliptic equations of interest here are the constraint 
equations that must be solved as part of the general 
relativistic Cauchy initial data problem. We focus on 
two axisymmetric problems: the initial data for a black 
hole spacetime with angular momentum, and a spacetime 
with a black hole superposed with gravitational waves 
(Brill waves) . Solutions for both of these problems have 
been found by others using FD : our focus here is to 
demonstrate the use of PSC for these problems and com- 
pare the computational cost of high accuracy solutions 
obtained using both PSC and FD. 

In section O we review briefly the key constraint equa- 
tions that arise in the traditional space-plus-time decom- 
position of the Einstein field equations and describe three 
different elliptic problems — a nonlinear model problem 
whose analytic solution is known, the nonlinear Hamilto- 
nian constraint equation for an axisymmetric black hole 
spacetime with angular momentum, and the Hamiltonian 
constraint equation for a spacetime with a black hole su- 
perposed with Brill waves — that have been solved using 
FD and that we solve here using PSC. We describe PSC 
in section and compare the computational cost of 
PSC and FD for high accuracy solutions in section |M. 
In section |^ we solve the problems described in section^ 
using PSC and compare the performance of PSC with 
that obtained by other authors using FD. In section [v|, 
we discuss our conclusions and their implications for solv- 
ing problems in numerical relativity. Finally, whether by 
FD or PSC the solution of a nonlinear elliptic system in- 
volves solving a potentially large system of (nonlinear) 
algebraic equations. We describe the methods we use for 
solving them in appendix 



II. INITIAL VALUE EQUATIONS 
A. Introduction 

The general relativistic Cauchy initial value problem 
requires that we specify the metric and extrinsic curva- 
ture of a three-dimensional spacelike hypersurface. These 
quantities cannot be specified arbitrarily: rather they 
must satisfy a set of constraint equations, which are a 
subset of the Einstein field equations. Arnowitt, Deser, 
and Misner (ADM) ||] were the first to formulate the 
Cauchy initial value problem in relativity in this way; 
however, the most common expression of this 3-1-1 de- 
composition is due to York 

In the York formulation of the ADM equations the 
spacelike hypersurfaces are taken to be level surfaces of 



some spacetime scalar function r. The generalized co- 
ordinates and conjugate momentum are the three-metric 
7ij induced on the spacelike hypersurface by the space- 
time metric and the extrinsic curvature Kij of the space- 
like hypersurface]^ Position on each surface is described 
by a set of spatial coordinates Xi (where Latin indices 
run from 1 to 3 and indicate spatial coordinates on the 
slice), so that the line-element on the hypersurface is 



(2.1) 



The normal n to the spacelike hypersurface is every- 
where timelike. The time coordinate direction t, how- 
ever, need not be exactly along the normal. We can write 
t in terms of n and the spatial vectors that span the tan- 
gent space of the hypersurface: 



d 



(2.2) 



where the lapse a describes how rapidly time elapses as 
one moves along the hypersurface normal n, and the shift 



^ ^ dx^ ' 



(2.3) 



is a vector field confined entirely to the hypersurface tan- 
gent space that describes how the spatial coordinates are 
shifted, relative to n, as one moves from one hypersur- 
face to the next. The lapse and shift are free functions: 
they correspond to the freedom to specify the evolution 
of the coordinate system that labels points in spacetime. 

The space-time line element at any point on a hyper- 
surface is related to the spatial metric at that point, the 
lapse, and the shift: 

ds^ ~ g^^dx'^dx'^ 

= - (a^ - PaP") dt^ + 2(3idx'dt + j.jdx'dx^ (2.4) 

(where Greek indices run from to 3 and include the 
time coordinate t, which is sometimes referred to as Xq). 

Given a spacetime foliation, choice of "field variables" 
"fij and Kij , and coordinates (embodied in the lapse and 
shift), the field equations can be decomposed into four 
constraint equations, which jij and Kij must satisfy on 
each slice, and six evolution equations, which describe 



*This approach is by no means unique. For example, recent 
work by Choquet-Bruhat, York and collaborators (cf. |lO|] and 
references therein) on a different choice of generalized coor- 
dinates and momenta have yielded approaches in which the 
evolution equations form a first-order symmetric and hyper- 
bolic (FOSH) system. Many powerful numerical methods for 
solving FOSH systems exist and the numerical relativity com- 
munity is only now beginning to explore how these solution 
techniques can be brought to bear on the field equations in 
this form. 
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how the three-metric and extrinsic curvature evolve from 
one shce to the next. In this paper we consider only the 
problem of consistent specification of initial data: i.e., 
the solution of the constraint equations. 

The four constraint equations (in vacuum) are 



where ^^■'Vi i s the covariant derivative associated with 
7ij . Equation 2.10a is generally referred to as the Hamil- 



^^^R + K"^ - KabK"'' = 0, 
(3)Va (K''' ~ if 7") = 0, 



(2.5a) 
(2.5b) 



where ('^^i? is the Ricci scalar associated with 7^ , '■'^•'Va 
is the covariant derivative associated with 7^ , and 



K := Kab^" 



(2.6) 



is the trace of the extrinsic curvature Ka. Note that, as 



befits constraints, equations 2_^ involve only derivatives 
of 7ij and Kij in the tangent space of the corresponding 
hypersurface. 



B. Conformal imaging formalism 

Our goal is to determine a 7^ and a Kij that satisfy the 
constraint equations and boundary conditions. These are 
both symmetric tensors on a three-dimensional hyper- 
surface; consequently, between the two there are twelve 
functions that must be specified. Equations 2.5 place 
only four constraints on these twelve functions. In order 
to solve the initial value problem, the remaining compo- 



tonian constraint, while equations 2.10b are generally re- 
ferred to as the momentum constraints. Since the only 
derivatives of Aij appear in the combination of a diver- 
gence, only the longitudinal part of Aij is constrained by 
the momentum constraints. 

Now turn to the boundary conditions. We are inter- 
ested in problems with a single black hole. Let the initial 
hypersurface be asymptotically flat, so that on the hy- 
persurface far from the black hole the curvature vanishes. 
Describe the black hole by an Einstein-Rosen bridge {i.e., 
by two asymptotically flat three-surfaces connected by a 
throat) and insist that the spacetime be inversion sym- 
metric through the throat. These choices impose the 
boundary conditions 



dip ip 
dr 2a 



lim ip{r) = 1 asymptotic flatness, (2.11a) 
= inversion symmetry, (2.11b) 



on tp where r = a is the coordinate location of the throat. 

We can now describe how to solve the constraint equa- 
tions. Let K vanish on the initial hypersurface; then the 
Hamiltonian and the momentum constraints (eqs. 2.10| ) 
decouple^ Pick a conformal background metric jij 
(which determines ^^^V) and transverse-traceless part 
Al^rp of the conformal extrinsic curvature. Solve the mo- 



nents of 7^ and K,j must be specified. York § has de- mentum constramts^(eqs 
veloped a convenient formalism, referred to as conformal 
imaging, for dividing the spatial metric and extrinsic cur- 
vature into constrained and unconstrained parts, which 
we summarize in this subsection. 

Associate 7^ with a conformal background three-metric 
jij through a conformal factor ip: 



(2.7) 



The extrinsic curvature K'^^ is split into its trace K and 
its trace-free part 



(2.8) 



The trace K is treated as a given scalar function which 
will be specified. The trace-free extrinsic curvature A^^ 
of the conformal metric 7^ can be expressed in terms of 
■0 and A'i-. 



A'^ =i/;^M'J. 



(2.9) 



2.10b ) for the longitudinal part 
of the trace-free conformal extrinsic curvature Aij. To- 
gether with Aj^rp the trace-free conformal extrinsic cur- 
vature is th us full y determined and the Hamiltonian con- 
straint (eq. 2.10a) can be solved for the conformal factor. 
This determines the three-metric 7^ and its extrinsic cur- 
vature K"^^ and completes the specification of the initial 
data. 



C. Three test problems 

1. Black hole with angular momentum 

Focus first on the initial data corresponding to an ax- 
isymmetric black hole spacetime with angular momen- 
tum. This problem was first examined analytically by 
Iprf , and has been explored numerically by [^|| . 

Choosing the conformal background metric to be flat 
{i.e., jij = Sij), [Tl| ] found an analytic solution to the 
momentum constraints (eqs. 2.10b) that carries angular 



The constraint equations (eqs. p.5| ) can also be ex- 
pressed in terms of the conformal background metric and 
its trace-free extrinsic curvature: 

- - + lAabA^'ip-' = 0, (2.10a) 

VaA''' - |v'7"Vaif = 0, (2.10b) 



momentum and obeys the isometry condition at the black 
hole throat. Draw through any point a sphere centered on 



^Such a hypersurface is said to have vanishing mean curva- 
ture, or be maximally embedded. 
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the black hole throat, and let be the outward-pointing 
unit vector to the sphere there. Letting J' be the an- 
gular momentum of the physical space, the Bowen-York 
solution to the momentum constraints is 



(2.12) 



Corresponding to this solution is the Hamiltonian con- 
straint (eq. 2.10a) for the conformal factor -0, 



^2 , 9 J2 sin^ e , 7 



0, 



together with its boundary conditions, 
lim ip{r) — 1, 

9-0 







r— >oo 

1 

dr ' Ya 

d'4' 



-0, 



^ -0, 



(2.13) 

(2.14a) 
(2.14b) 

(2.14c) 



= 0, TT 



which result from asymptotic flatness, the isometry con- 
dition at the throat (which is located at r = a), and 
axisymmetry respectively. The equation with boundary 
conditions for the conformal factor is solved on the do- 
main a < r < oo. Once the conformal factor is de- 
termined, the geometry of the initial slice is completely 
specified. 

A useful diagnostic of an initial data slice is to compute 
the total energy contained in the slice. O Murchadha 
and York |lj] have examined the ADM energy (cf . ||] ) in 
terms of the conformal decomposition formalism giving 



Eadm = Eadm — 



1 

2^ 



(2.15) 



where Eadm is the energy of the conformal metric. Thus 
when the conformal metric is flat, the ADM energy re- 
duces to 



1 

2^ 



Eadm = -— <t> ^'i^d^S,, 



(2.16) 



i.e., it is proportional to the integral of the normal com- 
ponent of the gradient of the conformal factor about the 
sphere at infinity. 



2. A model problem 

Bowen and York |jll[ also describe a nonlinear "model" 
of the Hamiltonian constraint equation that can be solved 
exactly, which we utilize in section ^ to test our code. 
The model equation is 



ip-' = 0, 



(2.17) 



with P a constant. Togeth er wi th the bou ndar y condi- 
tions described above (eqs. 2.14 ), equation 2.17 has the 
solution 



where 



2E 

— + 6- 

r 



E = 



2a^E a" 



41 1/4 



4a 



2\l/2 



(2.18a) 



(2.18b) 



If we evaluate equation 2.16| for this solution, we find that 
it has ADM energy E. 



3. Black hole plus Brill wave 

The second physical problem upon which we demon- 
strate the use of spectral methods for numerical relativity 
is that of a black hole superposed with a Brill |l^ wave, 
a problem studied using FD by 0. 

Following 0, let the initial slice be a spacetime isom- 
etry surface (i.e., time symmetric); then, the extrinsic 
curvature Ka vanishes and the momentum constraints 



(eqs. 2.10b) are trivially satisfied. 

To determine the conformal factor the Hamiltonian 



constraint (eq. 2.10a) must be solved, which requires the 
specification of a conformal background metric. Let the 
line-element of the conformal background metric have the 
form 



ds' 
where 



= [e^i {dr 



q := Asm" 



exp 



V In-, 
a 



^d9^ 



exp 



9 • 9 

r sm 



(2.19) 



G 



(2.20a) 
(2.20b) 



n is an even integer, and A, ryo, and tr are constant pa- 
rameters that describe the superposed Brill wave's ampli- 
tude, position, and width, respectively. With this choice 
the Hamiltonian constraint equation becomes 



(920 2 50 



1 cot 61 30 



dQ 



r-2 9612 
d^q 1 
9r2 r dr 



1 



(2.21) 



The boundary conditions of asymptotic flatness and 



inversion symmetry are again given by equations 2.14 



Furthermore, since the conformal metric has no "1/r" 
parts in its expansion at infinity, its energy vanishe s and 
the ADM energy for these slices is given by equation 2.16. 
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III. SPECTRAL METHODS 

A. Introduction 

Consider an elliptic differential equation, specified 
by the operator L on the d-dimcnsional open, simply- 
connected domain 2?, with boundary conditions given by 
the operator S on the boundary dV: 



L{u){x) = f{x) 
S{u){x)^g{x) 



xev, 
X e dv. 



(3.1a) 
(3.1b) 



There may be more than one boundary condition, in 
which case we can index S and g over the boundary con- 
ditions. 

Approximate the solution u(x) to this system as a sum 
over a sequence of basis functions (t>k{x) on 2? + dV, 



un{x) 



N-l 

^ Uk(j)k{x), 
k=0 



(3.2) 



where the Uk are constant coefficients. Corresponding to 
the approximate solution ujv is a residual Rn on 2? and 
rjv on dV: 

Rn = L{um) — f on 2?, (3.3a) 
rjv = S{un) — g on dV. (3.3b) 

The residual vanishes everywhere for the exact solution 
u. 

In PSC we determine the coefficients Uk by requiring 
that wjv satisfies the differential equation and boundary 
conditions exactly at a fixed set of collocation points a;„: 
i.e., we require that 



= L[uN{Xn)] - fiXn) 
= S[uN{Xn)] - giXn) 



for Xn in 2?, 
for Xn on 92?, 



(3.4a) 
(3.4b) 



for all n. When the expansion functions and colloca- 
tion points are chosen appropriately a numerical solution 
of these equations can be found very efficiently. In the 
following subsection we discuss choices of the expansion 
basis and collocation points. 



B. Expansion basis and collocation points 

In PSC we require that the approximate solution upf 
satisfies the differential equation and boundary condi- 
tions exactly at the N collocation points a;„. The basis 
(/)/c should not constrain the values of the approximation 
at the collocation points; correspondingly, we can write 
the basis as a set of N functions (j)k{x) that satisfy a dis- 
crete orthogonality relationship on the collocation points 



N-l 



(t>jiXnWkiXn) 



(3.5) 



where the Uk are normalization constants. Note that the 
basis functions are inextricably linked with the colloca- 
tion points. 

It is sometimes the case that the basis can be chosen so 
that the boundary conditions are automatically satisfied. 
For example, consider a one-dimensional problem on the 
interval 



I =["1,1]. 



(3.6) 



If the boundary conditions are periodic then each element 
of the basis 



(j)kix) = exp [tti {x -I- 1) fc] , 



(3.7a) 



satisfies the boundary conditions; correspondingly, the 
approximate solution un automatically satisfies the 
boundary conditions. If, in addition, we choose the col- 
location points 



_ 2n 



(3.7b) 

then the basis satisfies the discrete orthogonality relation 



Sjk = — 



1 



N ^ 

n=0 



<t>j{Xn)(j>k{Xn 



(3.7c) 



In an arbitrary basis, or with arbitrarily chosen collo- 
cation points, finding the Uk from the UN{xn) requires 
the solution of a general linear system of N equations 
in N unknowns, which involves 0{N^) operations. For 
the basis and collocation points given in equations 3/7 
the ilk can be determined from the UN{xn) quickly and 
efficiently via the Fast Fourier Transform in 0{NhiN) 
operations. 

Arbitrary derivatives of the uat can also be computed 
quickly: writing 



(Pun 
dxP 



N-l 



fc=0 



we see immediately that 



{TTikf 



(3.8a) 



(3.8b) 



n=0 



Consequently, any derivative of ujv can be evaluated at 
all the collocation points in just 0{N In N) operations. 

The ability to evaluate efficiently the derivatives of 
UN at the collocation points is much more important 
than finding a basis whose individual members satisfy the 
boundary conditions. In the case of periodic boundary 
conditions we can have our cake and eat it, too. More 
generally we choose a basis in which we can efficiently 
compute the derivatives of mat at the collocation points 
and require separately that the approximate solution u n 
satisfy the boundary conditions at collocation points on 
the boundary. 
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For general boundary conditions a basis of Chebyshev 
polynomials often meets all of our requirements. ^ Recall 
that the Chebyshev polynomials are defined on I by 



Tk{x) = cos (fc cos ^ x) 



(3.9) 



A simple recursion relation allows us to find the deriva- 
tive of UN as another sum over Chebyshev polynomials: 
if I 



N 



UN 



then 



duN 
dx 



fe=0 



N~l 



(x) = u'tJkix), 



where 



with 



CfcMfc = Ufc^_2 + 2(fc + l)uk+l 



Ck 



2 fc 
1 fc > 1. 



(3.10) 



(3.11) 



(3.12) 



(3.13) 



If we choose collocation points x„ (for < n < N) 
according to 



N 



(3.14) 



then the Chebyshev polynomials satisfy the discrete or- 
thogonality relation 



Sjk 



Nck c 



—Tj{Xn)Tk{Xn), 



where 



Ck 



2 fc = or TV 
1 otherwise. 



(3.15) 



(3.16) 



Finally, exploiting the relation between the Chebyshev 
polynomials and the Fourier basis (cf. 3.9) allows us to 
find the Uk from the UN{xn) in 0{N\nN) time using a 



*The geometry of a problem might suggest other expan- 
sion functions, such as Legendre polynomials; however, a 
Chebyshev expansion does quite well and has the added con- 
venience that, with appropriately chosen collocation points, 
only C'(A''lnA'') are required to convert from the expansion 
coefficients to the function values at the collocation points 
and vice versa 

^For Chebyshev bases the conventional notation is that k 
runs from to A'^, not — 1; thus, there are -I- 1 coefficients 
and collocation points. 



fast transform appendix B]. With an expansion ba- 
sis of Chebyshev polynomials and an appropriate choice 
of collocation points we can thus evaluate derivatives of 
arbitrary order at the collocation points in 0{NhiN) 
operations. 

For problems on an arbitrary domain of dimension d 
greater than unity it is rarely the case that we can find 
a basis which permits rapid evaluation of derivatives. If 
the domain can be mapped smoothly to I"^ then we can 
write 

Ar(i) jv*'') 

UNir>...Nw{x) ^ X! "' X! '^k^■■■kAk^■■■kA^)^ (3.17a) 



where 



(3.17b) 
(3.17c) 



and the 



kit) 



for fixed is a basis on I which permits 



fast evaluation of derivatives with respect to its argument 
{e.g., Chebyshev polynomials). Associated with each set 
of basis functions are the collocation points xh. ; corre- 
spondingly, the collocation points associated with <f>ki---kcL 
are just the Ni - ■ ■ A^^-tuples 



(4V 



(3.17d) 



With this choice of basis and collocation points we can 
evaluate efficiently arbitrary derivatives of an approxi- 
mation Unm...nw ■ If the domain cannot be mapped 
smoothly to a c?-cube in , either more sophisticated 
methods such as domain decomposition must be 

used, or the problem may not be amenable to solution 
by PSC. 



C. Solving the system of equations 

The expansion basis, collocation points and differen- 
tial equation with boundary conditions determine a sys- 
tem of equations for the coefficients Uk or, equivalently, 
the approximate solution u n evaluated at the collocation 
points. Iterative solution methods (which require as few 
as C(A^lniV) operations) work well to solve the kind of 
systems of equations that arise from the application of a 
PSC method. 

If the elliptic system being solved is linear then the 
algebraic equations arising from either a FD or a PSC 
method are also linear and a unique solution is guaran- 
teed. If, on the other hand, the differential system is 
nonlinear, then the equations arising from FD or PSC 
are also nonlinear and a unique solution is not guaran- 
teed. Newton's method ||l|, sec. 12.13 and appendices C 
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and D] , where one solves the linearized equations begin- 
ning with a guess and then iterating, works well for these 
types of equations. As long as a good initial guess is 
chosen, the iteration will usually converge. In appendix 
^ we describe in detail the variant of Newton's method 
(Richardson's iteration) that we have used to solve the 
nonlinear system of algebraic equations that arise when 
we apply PSC to solve the Hamiltonian constraint equa- 
tions as posed in section 



IV. COMPARING FINITE DIFFERENCE AND 
PSEUDOSPECTRAL COLLOCATION METHODS 

A. Introduction 

Finite differencing and pseudospectral collocation are 
alternative ways to find approximate solutions to a sys- 
tem of differential equations. Consider the Poisson prob- 
lem in one dimension: 



The FD solution to equations 4.1 begins by approxi 



fix), (4.1a) 

on the interval I with Dirichlct boundary conditions 

w(-l) = u{l) = 0. (4.1b) 

In a FD approach to this problem we seek the values of 
u at discrete points a;„, say 



_ 2n 



for n — 0,1, . . . N . Algebraic equations are found by ap- 
proximating the differential operator cPu/dx'^ in equation 
4.1a by a ratio of differences: e.g., 



d u 



Un+l - 2u„ + Un- 

Ax^ 



for integer rt = 1, 2, . . . — 1 where 
Ax := 2/N. 



(4.3) 



(4.4a) 
(4.4b) 



4.1a 



With this discretization the differential equation 
yields — 1 equations f or the A^ + 1 unknown ii„. The 
boundary conditions (eq. 4.1b) yield two more equations, 
completely determining the it„: 



Un+l - 2Un + Un- 

Ax^ 



for 



1,2,. ..N-1, 
(4.5a) 



mating the differential equations. In the PSC method, 
on the other hand, we first approximate the solution at 
all points in I by a sum over a finite set of basis functions. 
For this example, we choose a Chebyshev basis; so, we 
write 



UNix) 



N 

E 

k=0 



UkTk{x). 



(4.6) 



Now insist that ujv satisfies the differential equation 
and boundary conditions exactly at the collocation points 



cos ■ 



N 



(4.7) 



for 71 = 0, 1, . . . A^. In particular, we require that the 
boundary conditions are satisfied and that, in addition, 
the differential equation is satisfied for integer n ranging 
from 1 to A^ — 1 : 



un{xo) 
unIxn) 



0, 
0, 



d^UN 

dx'^ 



{Xn) = f{Xn)- 



(4.8a) 
(4.8b) 

(4.8c) 



To evaluate equation 4.8c note that d'^u^/dx'^ can be 
written as 



d'^UN 

dx'^ 



(Xn) 



N 

E 

m=0 



dl^luNiXm). 



(4.2) The dnm can be determined by noting that 



d'^UN 

dx^ 



with 



and 



Cku'k 



(Xn) 



^+2 



N 

J2ukTk{Xn), 

k=0 



2{k + l)uk+i, 



Uk 



2 



^ 1 
Si 

n=0 



UN{Xn)Tk{Xn), 



(4.9a) 



(4.9b) 



(4.9c) 



(4.9d) 



where Ck and Ck are given by equations 3.13 and 3.16, 
respectively. 

The result is, again, a set of algebraic equations for 
MAr(x„): the values of the approximate solution at the 
collocation points. Finding the un{x,i) yields an approx- 
imate solution to the differential equation over the entire 
domain I sinc e the spectral coefhcients Uk are given by 
equation 



Mo = 0, 

UN — 0. 



(4.5b) 
(4.5c) 



The solution to these equations is the FD approximation 
to u{x) at the points Xn- 



**Alternatively, we could have constructed a system of equa- 
tions in terms of the unknown spectral coefficients. This 
would correspond to a spectral tau method: cf. Mt 
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For the linear problem posed here the solution to the 
algebraic system of equations that arise in either a FD 
or PSC solution can be solved directly or by any of the 
many standard iterative methods. For nonlinear prob- 
lems the systems are generally solved by linearizing the 
equations about an initial guess and then iterating the 
solution until it converges. We discuss one method of 
solution in appendix 



B. Convergence of approximations 

In either a FD or PSC solution to a differential equa- 
tion with boundary conditions we expect that, as N tends 
to infinity, the approximate solution should become ar- 
bitrarily accurate. For large N, the L2 error in a FD ap- 
proximation converges upon the exact solution as N^p 
for positive integer p. The value of p depends on the 
smoothness of / and the error in the approximation of 
the differential operator (in the example above, /dx^). 
Assuming that / is smooth the rate of convergence (mea- 
sured by the L2 error of the FD solution) is N~p when the 
truncation error of the differential operator is 0(Aa:^). 

In contrast, when the solution u is smooth the error 
made by a properly formulated spectral approximation 
decreases faster than any fixed power of N (where N 
is now the number of collocation points or basis func- 
tions) .0 For a heuristic understanding of this rapid 
convergence, note first that a PSC solution's deriva- 
tives at each collocation point involve all the {wAr(x„)} 
(cf. eq. 4.9). Correspondingly, it is as exact as possi- 
ble, given the information available at the N colloca- 
tion points. This suggests that an order N collocation 
spectral approximation to the derivatives of the unknown 
should make errors on order 0{/S.x^). The interval Ax, 
however, is also proportional to A^^^; so, we expect that 
the error in the spectral solution ujv should vary as 
0{N~^). A more rigorous analysis using convergence 
theory chapter 2] shows that for any function which 
is analytic on the domain of interest, a Chebyshev ex- 
pansion will converge exponentially (i.e. as 0{e^^)). If 
the function is also periodic then a Fourier expansion will 
converge exponentially. 



C. Computational cost of solutions 

The computational cost, in time, of a FD solution to 
a system of elliptic differential equations scales linearly 
with the number of grid points N while the accuracy e 
of the solution scales as A^~^', where p is the order of the 



^^In addition the individual spectral coefficient Uk should 
decrease exponentially with A'^ once the problem is sufficiently 
resolved. 



FD operator truncation error. Correspondingly, the cost 
-fCpD for a given accuracy scales as 



(4.10a) 



The cost i^psc of a PSC solution to the same system, on 
the other hand, scales as A^ In A^ (for an iterative solution) 
while e scales as exp(— A); consequently, the cost scales 
with accuracy e as 



ATpsc ^ -(Ine) In (- Ine) . 



(4.10b) 



Since it is the computational cost required to achieve a 
given accuracy that is important, the more rapid conver- 
gence of a PSC solution confers upon it a clear advantage. 
This advantage is made clear by considering how the ra- 
tio of costs scales with accuracy: 



ATpsc 



FD 



-e^/'P Ine In (-Ine) , 



(4.11) 



which tends to zero with e; consequently, increasing ac- 
curacy with a PSC solution is always more efficient than 
with a FD solution. 

The equations that arise from either a FD or PSC 
treatment of an elliptic differential system are typically 
solved using iterative methods; thus, at fixed resolution 
the storage requirements for either solution method are 
equivalent. As we have seen, however, fixed resolution 
does not correspond to fixed solution accuracy. As the 
desired solution accuracy increases the storage require- 
ments of a PSC solution fall relative to those of an FD 
solution a factor of — e^/^ Ine. 



V. SOLVING THE HAMILTONIAN 
CONSTRAINT 

A. Nonlinear model problem 

As a first example we solve the mo del Ha miltonian 
cons traint equation described in section II C 2 , equation 
ITtI: 



474 



r 



(5.1a) 



for r G [a, oo) with P a constant and with the boundary 
conditions 



or 2a 



86 



lim 7/'(r) = 1 asymptotic flatness, (5.1b) 
= inversion symmetry, (5.1c) 

= axisymmetry. (5. Id) 



?=0,7r 



For this model problem the solution ip is 



8 



where 



2E a 
— +6^ 

r r 



E = 



2a^E 



1/4 



(5.2a) 



(5.2b) 



is also the ADM energy for the initial data (cf. eq. ^.16 ). 

As described this problem is spherically symmetric; 
nevertheless, we treat it as axisymmetric to illustrate the 
methods used to solve the Hamiltonian constraint for the 



black hole with angular momentum (cf. sec. II C 1 ) a nd 
the black hole plus Brill wave problems (cf. sec. pC^ ). 

As a first step we map the domain r G [a, cxd), 9 G [0, tt] 
to a square in R^: letting 



2a 

X — 1, 

r 

y = cos 6, 



(5.3a) 
(5.3b) 



we have x G (—1, 1] and y G [—1, 1]. In terms of the (x, y ) 
coordinates, the model Hamiltonian constraint (eq. 5.1a) 
becomes 



subject to the boundary conditions 



lim ip = 1, 



9-0 
dx 



= 0. 



(5.5a) 
(5.5b) 



Note that with our choice of variables a nd ex pansion 
bases the angular boundary conditions (eq. 5. Id) are au- 
tomatically satisfied. 

Since ip is not periodic in either x or y, we adopt a 
Chebyshev basis for the approximate solution: 



j=0 k=0 

with the corresponding collocation points 

x,=cos— , 
Trfc 



Vk 



V 



For this problem, focus on approximations 



(5.6a) 

(5.6b) 
(5.6c) 

(5.7) 



for integer £. We keep Ny fixed as the model problem is 
independent of y. 



Following the discussion in appendix solve the PSC 
equations using Richardson's iteration with a second- 
order FD preconditioner. To obtain VP^, we need an initial 
guess '3/1°^ to begin the iteration. For the lowest resolu- 



tion expansion {N^ 
guess 

,(0) 



4) begin the iteration with the 



(3 



(5.8) 



which is the trivial solution for P = 0. Applying Richard- 
son's iteration will then give us the approximate solu- 
tion ^fi. Through the expansion 5.6a this determines an 
approximation for ip everywhere; in particular, it deter- 
mines an approximation at the collocation points corre- 
sponding to = 8, which we then use as the initial 
guess for determining the approximate solution In 
this same way we use a lower-resolution approximate so- 
lution as the initial guess for the approximate solution at 
the next higher- resolution, i.e. 



(5.9) 



To investigate the accuracy of our solution as a func- 
tion of resolution (basis dimension for PSC, number of 
grid points for FD) we evaluate a number of solutions 
differing only in resolution and evaluate several different 
error measures. 

1. For this problem we know the exact solution (cf. 
5.2); so, we calculate the L2 norm of the absolute 
error as a function of £: 




1 



j=0 k=0 



N^NyCjCk 



["itixj.Xk) 



1/2 



-tl}{Xj,Xk)Y 



(5.10) 



where Ck is given by equation 3.16 



2. We can also characterize the convergence of the ap- 
proximate solutions ^1 by calculating the L2 norm 
of the difference between the successive approxi- 
mate solutions: 



5^p 



I* 



2 • 



(5.11) 



The errors 5"^ and A^P are defined for either FD or 
PSC solutions. 

3. We also evaluate, by analogy with A^E* and (5^, the 
quantities 



AEi = \Ei~E\, 
SEg = \Ei — Ee-i\ 



(5.12) 
(5.13) 



where Ei is the ADM mass-energy associated with 
the approximate confo rmal factor '^g. We evaluate 
E using equation 2.16. 
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FIG. 1. Spectral convergence for a nonlinear model prob- 
lem. Plotted are a measure of the absolute error A^'^ , and two 
approximate measures of the error S'i'e and S'^'e as a function 
of Nx, the number of radial functions, for the case P = 1. 



FIG. 2. Spectral convergence for the solution of the Hamil- 
tonian constraint equation for a black hole with angular mo- 
mentum . Plotted are three approximate measures of the 
error S'ile, 5^t and SE as a function of Nx, the number of 
radial functions, for J = 1. 



4. For PSC solutions only we define the relative error 
measure 



j=0 k=0 



(5.14) 



with an accuracy of AE « 10^^*^ is obtained by doubling 
the number of radial functions. To achieve the same ac- 
curacy the FD approximation would require (assuming 
second order FD) a resolution of 3 x 10^ radial points. 



which characterizes the changes in the spectral co- 
efficients as the order of the approximation in- 
creases. 

For a properly formulated spectral method, all of our 
error measures should decrease exponentially with N if 
the solution to the problem is analytic. 

Figure |l| shows the absolute and relative errors A'i'e 
and S'i'i, along with the change in the spectral coeffi- 
cients S'^g, for P = 1. The exponential convergence of 
the solution with increasing is apparent. Experience 
shows that as the problem becomes more nonlinear (i.e., 
P becomes larger) more terms are needed in the expan- 
sion in order to achieve the same accuracy. 

This system of equations has also been solved using 
FD methods A point comparison is telling: in ||] 
a second order accurate FD solution with a resolution 
of 1024 radial points points were required for a solution 
with a. AE ~ 10-^ independent of P. The PSC solu- 
tion described here achieves the same accuracy using an 
expansion with only 12 radial functions for P = 1, and 
24 functions for P = 10. In either case a PSC solution 



B. Black hole with angular momentum 

Now turn to consider a truly non-radial, but still ax - 
isymmetric, problem: a rotating black hole (cf. II C 1 ). 
As before (cf. 5.3) we map the semi-infinite domain r > a 
to the finite box x G (—1,1], y £ [—1,1], obtaining the 
system of equations 



9 / J 



dy"^ dy 



+ 64 (0: + 1)^(1 -2/2)^-^ = 0, (5.15) 

sub ject to the boundary conditions given in equations 

For this problem we do not have the exact solution; 
so, we consider only the relative errors ^VP, 6"^ and SE. 
Figure || (||) shows these quantities as functions of Nx for 
J/M^ equal to 1 (100). For these solutions ^'^ = '04£,Ar , 
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FIG. 3. Same as figure g with J = 100. 

where initially Ny — 4 and is incremented by two^ when- 
ever the difference between S'^i with and without the 
increment was greater than ten percent. Again we see 
rapid, exponential convergence of the solution with . 

This problem has also been solved using second order 
FD 1^ . For a solution accuracy SE ~ 10-^ § found that 
a resolution 1024 radial and 384 angular grid points was 
required, roughly independent of the value of J. We find 
that PSC achieves the same accuracy with an expansion 
basis of 12 radial (and 4 angular) functions for J — 1, and 
24 radial (and 8 angular) functions for J ~ 100. Solution 
accuracies of 10^^" can be obtained for the PSC solution 
simply by doubling the size of the expansion basis (in 
X and y). For a similar increase in accuracy of the FD 
solution a grid approximately 300 times larger in each 
dimension would be required. 



C. Black hole plus Brill wave 

As a final example we consider the Hamiltonian con- 
straint for a black hole superposed with a Brill wave. Af- 
ter mapping this problem to the {x, y) domain we obtain 



"'■''■Along with axisymmetry, this problem has equatorial 
plane symmetry so '^i is even in y. By exploiting this sym- 
metry, we could reduce our number of angular functions by a 
factor of two 



FIG. 4. Spectral convergence for the solution of the Hamil- 
tonian constraint equation for a black hole plus Brill wave. 
Plotted is an approximate measure of the error 5^t as a 
function of N^, the number of radial functions, for the case 
yl = ryo = cr = 1, n = 2. 



the system of equations 



22/ 



dy 



= 0, 



with 



R = {x + iy 



dx"^ 



dx 



.2\ ^ 

dy'^ 



(5.16a) 



dq 
dy 
(5.16b) 



where q is given by equation 2.20a, and subject to the 
boundary conditions 5.5. 

In figure ^ we show S'^ as a function of N.^., for the 
Brill wave parameters a — A = rjo = I and n — 2. 
For these solutions "ifi — ij}u,N where initially Ny — 
4, and is incremented by two whenever the difference 
between 5'^i with and without the increment was greater 
than ten percent. The convergence, while rapid, is not 
quite exponential. In addition, the nearly exponentially 
decreasing error is impressed with a wave that is nearly 
periodic in spectral resolution log N^. . We attribute this 



behavior to the resolution of the factor R (cf. eq. 5.16b 
and also eq. 2.20a for q). Figure^ shows the the error Ai? 
obtained when we form approximate Rn^,n according 
to 
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60 80 



FI G. 5. The error in the spectral representation of R (equa- 



tion 5.16b) for the case shown in figure 



with 



RN^,Ny -EE RjkTj{x)T,{y), (5.17a) 

j=0 fc=0 



4 ^ 

X y h. J j^^p t m 

(5.17b) 



The structure in the sohition is the same as the structure 
in the Chebyshev approximation to R. 

This problem has also been solved using FD methods 
, enabling us to compare the resolution required for ap- 
proximate FD or PSC solution for a given accuracy. With 
second order FD a solution whose error SE is 3 x 10^''' 
required a resolution of 400 radial and 105 angular grid 
points. To achieve the same accuracy the PSC solution 
described here requires a basis of only 36 radial (and 12 
angular) Chebyshev polynomials. 



VI. DISCUSSION 

Pseudospectral collocation is a very efficient way of 
solving the nonlinear elliptic equations that arise in nu- 
merical relativity. These problems typically have smooth 
solutions; correspondingly, the approximate solutions ob- 
tained using pseudospectral collocation converge upon 



the exact solution exponentially with the number of col- 
location points. As a result, the cost of a high accuracy 
pseudospectral solution is not significantly greater than 
the cost of a similar solution of moderate accuracy. Since 
the computational burden of solving the pseudospectral 
collocation equations with a given number of collocation 
points is no greater than that required to solve the finite 
difference equations for the same number of grid points, 
the computational demands of a pseudospectral colloca- 
tion solution are far less than those of a finite difference 
solution for even modest accuracy. 

Another important advantage of a pseudospectral col- 
location solution over a finite differencing solution in- 
volves the formulation of the boundary conditions. In a 
finite difference solution the boundary conditions must 
be reformulated as finite difference equations or incor- 
porated approximately into the formulation of the fi- 
nite difference operator of the differential equations be- 
ing solved. This generally involves the introduction of 
auxiliary boundary conditions, which are not part of the 
original problem. For example, consider the second order 
elliptic equation on I: 



u(-l) = ,i(l) 



0. 



(6.1a) 
(6.1b) 



A fourth-order accurate finite difference approximation 
to the differential operator (P /dx"^ is 



16(wj-i - 2uj + Uj+i) - (uj-2 — 2uj + Uj+2) 



where 



Uj = u{jAx). 



(6.2) 
(6.3) 



Befo re t his finite difference operator can be used in equa- 
tion 3.1 it must be modified at the grid points —1-1- Ax 
and 1 — Ax since —1 — Ax and 1 -I- Ax both lie outside the 
computational domain. In this case, four boundary con- 
ditions are required (at x equal to —1, —1-1- Ax, 1 — Ax 
and 1) even though the second order equation 6.1a prop- 
erly admits of only two boundary conditions. 

In pseudospectral collocation, on the other hand, no 
auxiliary boundary conditions need be formulated: the 
approximate solution is expressed as an analytic func- 
tion, which is required to satisfy the boundary condition 
equations exactly at the collocation points on the bound- 
ary. 

These advantages of pseudospectral collocation solu- 
tion come at a cost. When properly implemented the 
computational expense of pseudospectral collocation may 
be considerably less than the expense of finite differenc- 
ing; however, the difficulty of implementation is greater. 
The efficient solution of the algebraic equations arising 
from pseudospectral collocation generally require the use 
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of sophisticated iterative methods. Additionally, the ex- 
act solution itself should be smooth on the computational 
domain in order that exponential convergence is attained. 
Finally, and perhaps most importantly, for problems of 
dimension greater than unity the computational domain 
must be sufRciently simple that it can be mapped to I'' 
or be decomposed into sub-domains that can each be 
mapped to I'' {e.g., an L-shaped region can be decom- 
posed into two regions, each of which can be mapped to 

12). 

We have not here investigated the application of pseu- 
dospectral collocation techniques to evolution problems. 
Pseudospectral collocation methods have been used to 
solve problems in other fields {e.g., fluid dynamics) with 
great success [^^. Our own experience in applying 
these techniques to evolution problems in numerical rel- 
ativity shows promise, but is not yet complete. 
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APPENDIX A: SOLVING THE 
PSEUDOSPECTRAL COLLOCATION 
EQUATIONS 

In this appendix we describe one method of solving 
the nonlinear equations that arise from applying a PSC 
method to a nonlinear elliptic system of equations 



L{u) 
S{u) 



f onV, 
g on dV. 



(Ala) 
(Alb) 



Choosing an expansion basis and corresponding col- 
location points, the PSC solution of these equations is 
fully characterized by the values oiu^ at the collocation 
points Xn'. from these the coefficients of the expansion 
and all the derivatives of the approximate solution can 
be determined. Write the values of the approximate so- 
lution um at the collocation points Xn as a vector U , 



Un = UN{Xn)- 



(A2) 



Corresponding to the approximate solution un is a resid- 
ual Rn on V and on dV: 



Rn = L{un) - / onV, 
tn = S{un) — g on dV. 



(A3a) 
(A3b) 



The residual vanishes everywhere for the exact solution 
u. Write the values of the residual at the collocation 
points Xn as a vector R, 



" [ r]s{xn) Xn on dV. 



(A4) 



The PSC solution U satisfies the algebraic equations 

R[U]^Q. (A5) 



Before describing how to solve equation A5 for a non- 
linear system {i.e., nonlinear L or S) we describe the 
method of solution for a linear system. 



When the system of differential equations Al is linear 
so is the system of algebraic equations A5. In this case 
we can write 



AU = F, 



(A6) 
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where A is a matrix and is a vector whose components 
take on the values of / and g evaluated at the collocation 
points in the domain V and its boundary dV. In PSC 
the matrix A is typically full. Direct solution methods 
require 0{N^) operations for such systems; for efficiency 
such systems are generally solved by iterative methods, 
which typically requires many fewer operations to find an 
accurate solution. 

A simple and effective iterative method for solving 
equation A€ is Richardson's iteration. Suppose we have 
a guess to the solution U of equation A6 . A better 
approximation to U is given by 

^ yW _^jjW^ (A7) 
where the residual ii'*' vector is given by 

and w is a relaxation parameter, which must be deter- 
mined. The optimal value of uj and the rate of conver- 
gence of the iterations depend upon the eigenvalues of A. 
For Richardson's iteration the optimal w is 



'opt 



A. 



(A9) 



where A^ax and Amin are the largest and smallest eigen- 
values of A. This choice minimizes the spectral radius 



An 



A. 



P = 



A. 



+ A„ 



(AlO) 



of the iteration matrix^G = / — loK. The convergence 
rate of the iteration is 



7^ = -Inp. 



(All) 



The reciprocal of TZ. measures the number of iterations 
required to reduce the error by a factor of e. 

Richardson's iteration is, by itself, not necessarily more 
efficient than a direct solution method. Consider, for 
example, the second-order differential equation 



u(-l) = u{l) = 



(-1,1), 



(Al2a) 
(A12b) 



(cf. also section fV). A PSC solution with a Chebyshev 
expansion basis leads to an operator A with a spectral 
condition number Amax/Amin that is 0{N'^). This gives 
a rate of convergence TZ ~ 0{N~'^); correspondingly, 
0{N'^) iterations are required to obtain a reasonable solu- 
tion. Since each iteration requires 0{N\nN) operations 
{i.e., it is asymptotically dominated by the cost of eval- 
uating the derivatives d?u/dx^ given the A'^ -I- 1 UN^Xn)) 
the total cost of obtaining a solution U is 0{N^\'a-N), 
which is slightly more expensive than a direct solution. 



We can speed the convergence of Richardson's iteration 
by solving an equivalent problem whose spectral condi- 
tion number is better behaved. Introduce the precondi- 
tioning matrix H and consider the equivalent system 



H-^AU = H^F. 



(A13) 



Now given an approximation V*-*-* to U, a better approx- 
imation is given by 



(AM) 



where i?*^*^ is given as before and uj' is related to the 
eigenvalues of the linear operator H~^A. 

In practice we never actually invert the preconditioning 
matrix H; instead we solve 



H 



-UJ 



(A15) 



for successive approximations. In order that this 
equation for successive approximations should converge 
rapidly we require a preconditioning matrix H such that 



equation A15 is inexpensive to solve, and 



• the spectral condition number k' of H ^A is close 
to unity. 

If is a good approximation of A~^ then the second 
condition will be satisfied; consequen tly, w e look for ap- 
is inexpensive 



proximations to A for which equation A15 
to solve. 

The operator A arises from a system of differential 
equations. For one-dimensional problems a low-order FD 
approximation to this operator (with grid points coinci- 
dent with the collocation points) gives rise to a banded 
system with a small number of bands close to the main 
diagonal. When this FD operator i s use d as the pre- 
conditioner the system of equations A15 can be solved 
efficiently using direct methods.^ 



For instance, in the example considered here (eq. Al2| ) 
we can set H to be the second-order accurate FD oper- 
ator corresponding to L. The eigenvalues of the precon- 
ditioned operator H^A are all in the range 1 < A^*^ < 
7r^/4: i.e., the spectral condition number is independent 
of A. In this case the optimal relaxation parameter is 



opt 



4 
7' 



(A16) 



and each iteration reduces the residual by a factor of 
approximately 7/3 (independent of N) [Q. The asymp- 
totic cost of finding a solution is thus proportional to 
the cost of evaluating the residual, 0{N In N), which is 



^^For more details on the use of FD operators as precondi- 
tioners for spectral problems see p^ . 
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much more rapid than solution via a direct method or 
Richardson's iteration without a preconditioner. 

For higher dimensional problems the FD precondi- 
tioner still leads to a banded system with a small num- 
ber of non-zero bands; however, some of those bands are 



can 



found far from the main diagonal and equation A15 
no longer be solved efficiently using direct methods. If N 
becomes so large that the cost of solving these equations 
with the FD preconditioner is too great, then the equa- 
tions for the successive approximations can themselves be 
solved iteratively, other preconditioners can be explored 
(cf. [l|,|^,[T6t), or the original equations can be solved us- 
ing another iterative technique, such as multigrid Q. For 
the problems considered in this paper N never became 
so large that a direct solution of equation A15 with the 
FD preconditioner wa s problematic. 

If the equations Al are nonlinear, the algebraic equa- 
tions satisfied by U are similarly nonlinear. Write the 
nonlinear equations as 



C{U) = F, 



(A17) 



where C is a nonlinear function of U. In order to solve 
this nonlinear system of equations, we apply Newton's 
iteration sec. 12.13 and appendices C and D]. For 



equation |A5| , Newton's iteration is 



(A18) 



where A^(i) is the linear operator that arises from lin- 
earizing A about V^*-* and i?*-*-* is the nonlinear residual 
given by 



fiW ^ c{U) - F. 



(A19) 



Equation A18 is a linear system to be solved at each step 
of Newton's iteration. In the same way as before we can 
introduce a preconditioner, in which case we have the 
nonlinear Richardson's iteration 



(A20) 



Here H is any suitable preconditioning matrix for A-j/(i) . 
For the problems solved in section ^ we used as a pre- 
conditioning matrix a second-order accurate FD operator 
corresponding to the derivative terms of C ignoring the 
nonlinear terms. (Equivalently we could have used the 
FD operator corresponding to the linearized operator, 
but for the problems we examined this was not neces- 
sary.) 

As a quick demonstration, consider the example prob- 
lem 



TABLE I. The values of the absolute error /S.upsc of 
a PSC calculation, as well as the absolute error of a sec- 
ond-order FD calculation Auf d for several values of A'^ for 
the example problem (equation A21). For A'^ > 16, the PSC 
solution is contaminated with roundoff errors. 



N 






4 


1.7 xlO"^ 


2.5 xlO"^ 


8 


3.2 xlO"* 


5.4 xlO"^ 


16 


6.9 xlO-i° 


1.3 xlO-2 


32 




3.3 xlO"^ 


64 




8.2 xlO"* 


128 




2.1 xlO~* 



PSC approximation on a Chebyshev basis. For this prob- 
lem we know the exact solution. 



u{x) = sin(7ra;). 



(A23) 



Table | lists Aupsc and Aupjj (cf. 5.10) for increasing 
A^ (number of grid points for the FD approximation, ba- 
sis dimension for the PSC approximation). The rapid 
convergence of PSC is apparent. The second-order FD 
solution requires 128 points to equal the moderate accu- 
racy of an eighth-order PSC solution. In order to match 
the high accuracy of the 16th-order PSC solution would 
require a second-order FD solution with 6.5 xlO"* points. 



(Pu 

—r -e"'[(l -7r2)sin(7r2;) -|-27rcos(7rx)] = 0, (A21) 
ax 

u{-l)^u{l)^0. (A22) 

We have evaluated approximate solutions to this problem 
using a second-order accurate FD approximation and a 
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